function ConsDist=IRFAdd_ConsDist(PP,SS,MODEL,EquJac,T,IRF)

%% Reconver the Consumption Deviation
HH          =   MODEL.EquBlock.HH;
DevInput    =   zeros(HH.Dim.In,T);
for ii=1:length(HH.Var.In)
    InVar           =   HH.Var.In{ii};
    if strcmp(InVar,'ValFunp')
        DevInput(HH.LocIdx.In.ValFunp,1:T-1)=   IRF.ValFun(:,2:T);
    else
        DevInput(HH.LocIdx.In.(InVar),:)    =   IRF.(InVar)(:,1:T);
    end
end
OutputDev   =   EquJac.Input*DevInput;


%% State and Distribution
%==========================================================================
% QQ
%==========================================================================
Policy_c_SS =   [SS.DistPolicy.V.c;SS.DistPolicy.Vhat.c];
Dev_Policy_c=   OutputDev(1:length(Policy_c_SS),:);
Policy_c    =   Policy_c_SS.*exp(Dev_Policy_c);
Policy_dc   =   Policy_c-Policy_c_SS;
%==========================================================================
% QW
%==========================================================================
TempQW_c_SS =   SS.TrProb.UnitTrMat.ReEval*SS.Dist_fzo.QW_fzo;
QW_c_SS     =   [(1-PP.XI)*TempQW_c_SS;PP.XI*TempQW_c_SS];
Dev_QW_c    =   OutputDev(1+length(Policy_c_SS):end,:);
QW_c        =   QW_c_SS + Dev_QW_c;

%==========================================================================
% Groups
%==========================================================================
f           =   repmat(SS.DistApp.fzo.State(:,1),[2,1]);
ID          =   repmat(SS.DistApp.fzo.State(:,2),[2,1]);
N_State     =   size(ID,1);
UnitNum     =   [SS.DistApp.fzo.UnitNum(1),...
                 PP.ExoState.IdioInc.N,PP.ExoState.Idio_FI.N,PP.ExoState.Idio_RI.N,...
                 2];
% Financial and Real Integration
ID_z        =   SS.FunApp.ExoState.State(ID,1);
ID_FI       =   SS.FunApp.ExoState.State(ID,2);
ID_RI       =   SS.FunApp.ExoState.State(ID,3);

FI          =   PP.ExoState.Idio_FI.Node(ID_FI);
RI          =   PP.ExoState.Idio_RI.Node(ID_RI);
% Net Wealth Quantiles
QuantList   =   (0.2:0.2:0.8)';%[0.25;0.50;0.75];
Quant_f     =   DistApp_Hist_Quantile(SS.DistApp.fzo.UnitNode{1},SS.Dist_fzo.Marginal_fzo.f,QuantList);
Quant_f_Num =   length(Quant_f)+1;

Mid_f       =   max(SS.DistApp.fzo.UnitNode{1}(cumsum(SS.Dist_fzo.Marginal_fzo.f)<0.5));
% Group Index
Ind         =   struct();

Ind.dom     =   FI==0;
Ind.ext     =   ~Ind.dom;

Ind.N       =   RI==0;
Ind.H       =   ~Ind.N;

Ind.rich    =   f>Mid_f;
Ind.poor    =   f<=Mid_f;

Ind.Quant_f =   cell(Quant_f_Num,1);
Ind.Quant_f_Label ...
            =   cell(Quant_f_Num,1);
for ii=1:Quant_f_Num
    switch ii
        case 1
            Ind.Quant_f{ii}     =   f<Quant_f(ii);
        case Quant_f_Num
            Ind.Quant_f{ii}     =   f>=Quant_f(ii-1);
        otherwise
            Ind.Quant_f{ii}     =   (f>=Quant_f(ii-1)) & (f<Quant_f(ii));
    end
    Ind.Quant_f_Label{ii} ...
            =   ['Q',num2str(ii)];
end

Ind.f_Fine  =   cell(SS.DistApp.fzo.UnitNum(1),1);
for ii=1:SS.DistApp.fzo.UnitNum(1)
    Ind.f_Fine{ii}  =   f==SS.DistApp.fzo.UnitNode{1}(ii);
end

Ind.z       =   cell(PP.ExoState.IdioInc.N,1);
for ii=1:PP.ExoState.IdioInc.N
    Ind.z{ii}   =   (ID_z==ii);
end

Ind.Group.NA_Fine=Ind.f_Fine;
Ind.Group.z =   Ind.z;
Ind.Group.RI=   {Ind.H,Ind.N};
Ind.Group.FI=   {Ind.ext,Ind.dom};
Ind.Group.NA=   {Ind.rich,Ind.poor};

TempList    =   {{'RI','FI'},{'RI','NA'},{'FI','NA'},...
                 {'RI','FI','NA'}};
for ii=1:length(TempList)
    TempGroup   =   TempList{ii};
    TempLabel   =   strjoin(TempGroup,'_');
    TempN       =   length(TempGroup);
    Ind.Group.(TempLabel) ...
                =   cell(1,2^TempN);
    Ind.Group.(TempLabel){1}...
                =   ones(size(Ind.H));
    for jj=1:TempN
        TempStart   =   2^(jj-1);
        TempInd     =   Ind.Group.(TempLabel)(1:TempStart);
        for nn=1:TempStart
            Ind.Group.(TempLabel){nn} ...
                        =   TempInd{nn} & Ind.Group.(TempGroup{jj}){1};
            Ind.Group.(TempLabel){TempStart+nn} ...
                        =   TempInd{nn} & Ind.Group.(TempGroup{jj}){2};
        end
    end
end
IndGroupList=   fieldnames(Ind.Group);
%% Different Statistics
TopSharePct =   [0.1;0.9;0.25;0.75;0.5];
[Gini_SS,TopShare_SS] ...
            =   DistApp_Hist_InequalityMeasure(Policy_c_SS,full(QW_c_SS),TopSharePct);
[CVar_SS,CMean_SS,~] ...
            =   DistApp_QW_QW2Moment(QW_c_SS,Policy_c_SS,[2]);
C_CV_SS     =   sqrt(CVar_SS)/CMean_SS;
[LogCVar_SS,LogCMean_SS,~] ...
            =   DistApp_QW_QW2Moment(QW_c_SS,log(Policy_c_SS),[2]);
LogC_Std_SS =   sqrt(LogCVar_SS);
[dC_Gini_SS,~] ...
            =   DistApp_Hist_InequalityMeasure(Policy_c_SS./Policy_c_SS,full(QW_c_SS),TopSharePct);
for tt=1:T
    % Dispersion: Std and Inter-Quantile Distance
    Var         =   DistApp_QW_QW2Moment(QW_c(:,tt),Dev_Policy_c(:,tt),[2]);
    Quant       =   DistApp_Hist_Quantile(Dev_Policy_c(:,tt),QW_c(:,tt),[0.10;0.90;0.25;0.75]);
    Disp_Unit   =   [sqrt(Var);Quant(2)-Quant(1);Quant(4)-Quant(3)];
    if tt==1
        Disp        =   zeros(length(Disp_Unit),T);
    end
    Disp(:,tt)  =   Disp_Unit;
    % Variance Decomposition
    for ii=1:length(IndGroupList)
        TempGroup   =   IndGroupList{ii};
        VarDec_Unit =   CrossSectionVarDec(Dev_Policy_c(:,tt),QW_c(:,tt),...
                                           Ind.Group.(TempGroup));
        if tt==1
            VarDec.(TempGroup)  =   zeros(length(VarDec_Unit),T);
        end
        VarDec.(TempGroup)(:,tt) ...
                    =   VarDec_Unit;
    end
    % Consumption over Wealth Level
    Temp_1      =   QW_c_SS.*Policy_c_SS;
    Temp_2      =   QW_c(:,tt).*Policy_c(:,tt);
    
    
    AvgDev_byNA_Unit ...
                =   zeros(Quant_f_Num,1);
    for ii=1:Quant_f_Num
        AvgDev_byNA_Unit(ii)    =   log( sum(Temp_2(Ind.Quant_f{ii}))/sum(QW_c(Ind.Quant_f{ii},tt)) ) - ...
                                    log( sum(Temp_1(Ind.Quant_f{ii}))/sum(QW_c_SS(Ind.Quant_f{ii})) );
    end
    if tt==1
        AvgDev_byNA     =   zeros(length(AvgDev_byNA_Unit),T);
    end
    AvgDev_byNA(:,tt) ...
                =   AvgDev_byNA_Unit;
    % Dispersion of Level: Gini and Top Shares
    [Gini_Unit,TopShare_Unit] ...
                =   DistApp_Hist_InequalityMeasure(Policy_c(:,tt),QW_c(:,tt),TopSharePct);
    [CVar_Unit,CMean_Unit,~] ...
                =   DistApp_QW_QW2Moment(QW_c(:,tt),Policy_c(:,tt),[2]);
    [LogCVar_Unit,LogCMean_Unit,~] ...
                =   DistApp_QW_QW2Moment(QW_c(:,tt),log(Policy_c(:,tt)),[2]);
    [dC_Gini_Unit,~] ...
                =   DistApp_Hist_InequalityMeasure(Policy_c(:,tt)./Policy_c_SS,QW_c(:,tt),TopSharePct);
    if tt==1
        Gini        =   zeros(length(Gini_Unit),T);
        TopShare    =   zeros(length(TopShare_Unit),T);
        C_CV        =   zeros(length(CVar_Unit),T);
        LogC_Std    =   zeros(length(LogCVar_Unit),T);
        dC_Gini     =   zeros(length(dC_Gini_Unit),T);
    end
    Gini(:,tt)      =   Gini_Unit-Gini_SS;
    TopShare(:,tt)  =   TopShare_Unit-TopShare_SS;
    C_CV(:,tt)      =   sqrt(CVar_Unit)/CMean_Unit-C_CV_SS;
    LogC_Std(:,tt)  =   sqrt(LogCVar_Unit)-LogC_Std_SS;
    dC_Gini(:,tt)   =   dC_Gini_Unit-dC_Gini_SS;
end

ConsDist    =   struct('Disp',Disp,'VarDec',VarDec,'AvgDev_byNA',AvgDev_byNA,'Gini',Gini,'TopShare',TopShare,...
                       'C_CV',C_CV,'LogC_Std',LogC_Std,'dC_Gini',dC_Gini,...
                       'SS',struct('Gini',Gini_SS,'TopShare',TopShare_SS,'TopSharePct',TopSharePct,'QQ_f',(1:Quant_f_Num)'));